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' Abstract 

Although the spike-trains in neural networks are mainly constrained by the neural dynamics itself, global 
temporal constraints (refractoriness, time precision, propagation delays, ..) are also to be taken into account. 
. These constraints are revisited in this paper in order to use them in event-based simulation paradigms. 

We first review these constraints, and discuss their consequences at the simulation level, showing how event- 
based simulation of time-constrained networks can be simplified in this context: the underlying data-structures are 
' strongly simplified, while event-based and clock-based mechanisms can be easily mixed. These ideas are applied 

to punctual conductance-based generalized integrate-and-fire neural networks simulation, while spike-response 
model simulations are also revisited within this framework. 

As an outcome, a fast minimal complementary alternative with respect to existing simulation event-based 
methods, with the possibility to simulate interesting neuron models is implemented and experimented. 
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^ 1 Introduction 



Let us consider the simulation of large-scale networks of neurons, in a context where the spiking nature of neurons 
activity is made explicit (Gerstner & Kistler, 2002), either from a biological point of view or for computer sim- 
ulation. From the detailed Hodgkin-Huxley model (Hodgkin & Huxley, 1952), (still considered as the reference 
but unfortunately intractable when considering neural maps), back to the simplest integrated and fire (IF) model, 
• ^ , a large family of continuous-time models have been produced, often compared with respect to their (i) biological 
plausibility and their (ii) simulation efficiency. 

As far as this contribution is concerned, we consider a weaker notion of biological plausibility: a simulation is 
biologically plausible if it verifies an explicit set of constraints observed in biology. More precisely, we are going 
to consider a few global time constraints and develop their consequences at the simulation level. It appears here 
that these biological temporal limits are very precious quantitative elements, allowing us on one hand to bound and 
estimate the coding capability of such systems and, on the other hand, to improve simulations. 



Simulation efficiency of neural network simulators 

Simulation efficiency is a twofold issue of precision and performances. See (Brette et al., 2007) for a recent review 
about both event-based and clock-based simulation methods. 

Regarding precision, event-based simulations, in which firing times are not regularly discretized but calculated 
event by event at the machine precision level, provides (in principle) an unbiased solution. On the reverse, it has 
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been shown that a regular clock-based discretization of continuous neural systems may introduce systematic errors, 
with drastic consequences at the numerical level, even when considering very small sampling times (Rudolph & 
Destexhe, 2007). 

Furthermore, the computational cost is in theory an order of magnitude better using event-based sampling 
methods (Brette, 2006), although this may be not always the case in practice (Morrison, Mehring, Geisel, Aerstsen, 
& Diesmann, 2005), as further discussed in this paper. 

However, using event-based simulation methods is quite tiresome: Models can be simulated if and only if the 
next spike-time can be explicitly computed in reasonable time. This is the case only for a subset of existing neuron 
models so that not all models can be used. An event-based simulation kernel is more complicated to use than a 
clock-based. Existing simulators are essentially clock-based. Some of them integrate event-based as a marginal 
tool or in mixtures with clock-based methods (Brette et al., 2007). According to this collective review, only the 
fully supported, scientifically validated, pure event-based simulators is MVASpike (Rochel & Martinez, 2003), the 
NEURON software proposing a well-defined event-based mechanism (Hines & Carnevale, 2004), while several 
other implementations (e.g.: DAMNED (Mouraud, Paugam-Moisy, & Puzenat, 2006), MONSTER (Rudolph & 
Destexhe, 2006)) exists but are not publicly supported. 

In other words, event-based simulation methods may save precision and computation time, but not the scientist 
time. 

The goal of this paper is to propose solutions to overcome these difficulties, in order to have an easy to use 
unbiased simulation method. 

Considering integrate and fire models. 

At the present state of the art, considering adaptive, bi-dimensional, non-Unear, integrate-and-fire model with con- 
ductance based synaptic interaction (as e.g. in (Destexhe, 1997; Brette & Gerstner, 2005; Rudolph & Destexhe, 
2007)), called "punctual conductance based generalized integrate and fire models" (gIF), presents several advan- 
tages: 

- They seem to provide an effective description of the neuronal activity allowing to reproduce several important 
neuronal regimes (Izhikevich, 2004), with a good adequacy with respect to biological data, especially in high- 
conductance states, typical of cortical in-vivo activity (Destexhe, Rudolph, & Pare, 2003). 

- They provide nevertheless a simplification of Hodgkin-Huxley models, useful both for a mathematical anal- 
ysis and numerical simulations (Gerstner & Kistler, 2002; Izhikevich, 2003). 

In addition, though these models have mainly considered one neuron dynamics, they are easy to extend to 
network structure, including synaptic plasticity (Markram, Liibke, Frotscher, & Sakmann, 1997; Pfister & Ger- 
stner, 2006). See, e.g. (Rauch, La Camera, Luscher, Senn, & Fusi, 2003) for further elements in the context of 
experimental frameworks and (Camera, Giugliano, Senn, & Fusi, 2008a, 2008b) for a review. 

After a spike, it is assumed in such integrate and fired models that an instantaneous reset of the membrane 
potential occurs. This is the case for all models except the Spike Response Model of (Gerstner & Kistler, 2002). 
From the information theoretic point of view, it is a temptation to relate this spurious property to the erroneous fact 
that the neuronal network information is not bounded. In the biological reality, time synchronization is indeed not 
instantaneous (action potential time-course, synaptic delays, refractoriness, ..). More than that, these biological 
temporal limits are very precious quantitative elements, allowing one to bound and estimate the coding capability 
of the system. 
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Taking time-constraints into account 



The output of a spiking neural network is a set of events, defined by their occurrence times, up to some precision: 

:F={---t^---},t\<tl<---<t^ <■■■ ,^i,n 
where if corresponds to the nth spike time of the neuron of index i, with related inter-spike intervals = 

In computational or biological contexts, not all sequences correspond to spike trains since they are con- 
strained by the neural dynamic, while temporal constraints are also to be taken into account (Cessac, H.Rostro- 
Gonzalez, Vasquez, & Vieville, 2008). This is the key point here: Spike-times are 
[CI] bounded by a refractory period r, 
[C2] defined up to some absolute precision 5t, while 

[C3] there is always a minimal delay dt for one spike to propagate to another unit, and there might be (depending 
on the model assumptions) a 

[C4] maximal inter-spike interval D such that either the neuron fires within a time delay < D or remains quiescent 
forever). For biological neurons, orders of magnitude are typically, in milliseconds: 
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The derivations of these numerical values are reviewed elsewhere (Cessac et al., 2008). This has several 
consequence. On one hand, this allows us ti derive an upper bound for the amount of information: 

TV log2 [ji) bits during T seconds, 
while taking the numerical values into account it means for large T, a straight-forward numerical derivation leads 
to about IKhits / neuron. 

On the other hand (Cessac, 2008; Cessac & Vieville, 2008), it appears that for generalized integrate and fire 
neuron models with conductance synapses and constant external currents, the raster plot is generically periodic, 
with arbitrary large periods, while there is a one-to-one correspondence between orbits and rasters. This last 
fact, and the fact that more general models such as Hodgkin-Huxley neuron assemblies can be simulated during a 
bounded period of time (Cessac et al., 2008) provides theoretical justifications for the present work. 

Wliat is the paper about 

We develop here the consequences of the reviewed time constraints at the simulation level. Section 2 shows how 
event-based simulation of time-constrained networks can be impacted and somehow improved in this context. 
Section 3 considers punctual conductance based generahzed integrate and fire neural network simulation, while 
section 4 revisits spike-response model simulation within this framework. These mechanisms are experimented in 
section 5, where the computer implementation choices are discussed. 

Since the content of this paper requires the integration of data from the hterature reused here, we have collected 
these elements in the appendix in order to ease the main text reading, while maintaining the self-completeness of 
the contribution. 

2 Event-based simulation of time-constrained networks. 

Clock-based and event-based simulations of neural networks make already use at different level of the global time- 
constraints reviewed here. See e.g. (Brette et al., 2007) for an introduction and a large review and (Rochel & 
Martinez, 2003; Brette, 2006; Rudolph & Destexhe, 2006; Morrison, Straube, Plesser, & Diesmann, 2007) for 
simulations with event-based or hybrid mechanisms. However, it appears that existing event-based simulation 
mechanisms gain at being revisited. 
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In event-based simulation the exact simulation of networks of units (e.g. neurons) and firing events (e.g. spikes) 
fits in the discrete event system framework (Rochel & Martinez, 2003) and is defined at the neural unit level by : 

-1- the calculation of the next event-time (spike firing), 
-2- the update of the unit when a new event occurs. 

At the network level, the following two-stage mechanism completely implements the simulation: 

-a- retrieve the next event-time and the related unit, 

-b- require the update of the state of this unit, inform efferent units that this unit has emitted an 
event, and update the related event-times, 

repeating -a- and -b- when ever events occur or when some bound is reached. 

This mechanism may also take external events into account (i.e., not produced by the network units, but by 
external mechanisms). 

Such a strategy is thus based on two key features: 

• The calculation of the next event-time, conditioned to the present state and to the fact that, by definition, no 
event is received in the meanwhile for each unit; 

• The "future" event-time list, often named "priority queue", where times are sorted and in which event-times 

are retrieved and updated. 

The goal of this section is to revisit these two features considering [C3] and [C4]. 

Let us consider in this section a network with N units, an average of C cormections per units, with an average 
number M < N of firing units. 

2.1 Event-time queue for time-constrained networks 

Although, in the general case, spike-times must be sorted among pending events, yielding a 0(log(M)) complexity 
for each insertion, there exists data structure allowing to perform retrieve/update operations in constant "0(1)" 
time. Several efficient data structures and algorithms have been proposed to handle such event scheduling task. 
They are usually based on heap-like structures (Rochel & Martinez, 2003) or sets of buffers associated with some 
time intervals (such as the calendar queue in (Brette, 2006)). Ring buffers with fixed time step are used in (Morrison 
etal., 2005). 

The basic idea of these structures is to introduce buckets containing event-times in a given time interval. In- 
dexing these buckets allows one to access to the related times without considering what is outside the given time 
interval. However, depending on the fixed or adaptive bucket time intervals, bucket indexing mechanisms and 
times list managements inside a bucket, retrieve/update performances can highly vary. 

Let us now consider [C3] and assume that the bucket time-interval is lower than dt, the propagation delay, 
lower than the refractory period r and the spike time precision St. If an event in this bucket occurs, there is at 
least a dt delay before influencing any other event. Since other events in this bucket occurs in a dt interval they 
are going to occur before being influenced by another event. As a consequence, they do not influence each others. 
They thus can be treated independently. This means that, in such a bucket, events can be taken into account in any 
order (assuming that for a given neuron the synaptic effect of incoming spikes can be treated in any order within a 
dt window, since they are considered as synchronous at this time-scale). 

It thus appears a simple efficient solution to consider a "time histogram" with dt large buckets, as used in 
(Morrison et al., 2005) under the name of "ring buffer". This optimization is also available in (Rochel & Martinez, 
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2003) as an option, while (Brette, 2006) uses a standard calendar queue, thus more general, but a-priory less tuned 
to such simulation. Several simulation methods take into account [C3] (e.g. (Lee & Farhat, 2001; Connolly, 
Marian, & Reilly, 2004)). 

The drawback of this idea could be that the buffer size might be huge. Let us now consider [C4], i.e. the 
fact that relative event-time are either infinite (thus not stored in the time queue) or bounded. In this case with 
D = lO^'^Ws (considering fire rate down to O.lHz) and 10~[^'^1 (considering gap junctions) the buffer size 
S = D/dt = 10"''^^, which is easily affordable with computer memories. In other words, thanks to biological 
order of magnitudes reviewed previously, the histogram mechanism appears to be feasible. 

If [C4] does not hold, the data- structure can be easily adapted using a two scale mechanism. A value of D, 
such that almost all relative event-time are lower than D is to be chosen. Almost all event-times are stored in the 
initial data structure, whereas larger event-times are stored in a transient calendar queue before being reintroduced 
in the initial data structure. This add-on allows one to easily get rid of [C4] if needed, and still makes profit of the 
initial data structure for almost all events. In other words, this idea corresponds to considering a sliding window 
of width D to manage efficiently events in a near future. This is not implemented here, since models considered in 
the sequel verify [C4]. 

The fact that we use such a time -histogram and treat the events in a bucket in any order allows us to drastically 
simplify and speed-up the simulation. 

Considering the software implementation evaluated in the experimental section, we have observed the fol- 
lowing. Event retrieval requires less than 5 machine operations and event update less than 10, including the 
on-line detection of [C3] or [C4] violation. The simulation kernel^ is minimal (a 10Kb C-n- source code), using 
a 0{D/dt + N) buffer size and about 0{1 + C + e/dt) ~ 10 — 50 operations/spike, thus we a small overhead 
e ^ 1, corresponding to the histogram scan. In other words, the mechanism is "0(1)" as for others simulation 
methods. Moreover, the time constant is minimal in this case. Here, we save computation time, paying the price in 
memory. 

Remarks 

The fact that we use such a time-histogram does not mean that we have discretized the event-times. The 
approximation only concerns the way how events are processed. Each event time is stored and retrieved with its 
full numerical precision. Although [C2] limits the validity of this numerical value, it is indeed important to avoid 
any additional cumulative round-off. This is crucial in particular to avoid artificial synchrony (Rochel & Martinez, 
2003; Morrison et al., 2005). 

Using [C3] is not only a "trick". It changes the kind of network dynamics that can be simulated. For instance, 
consider a very standard integrate and fire neuron model. It can not be simulated in such a network, since it can 
instantaneously fire after receiving a spike, whereas in this framework adding an additional delay is required. Fur- 
thermore, avalanche phenomena (the fact that neurons instantaneously fire after receiving a spike, instantaneously 
driven other neurons and so on..) cannot occur. A step further, temporal paradoxes (the fact, e.g., that a inhibitory 
neuron instantaneously fires inhibiting itself thus . . is not supposed to fire) cannot occur and have not to be taken 
into account. When considering the simulation of biological systems, [C3] indeed holds. 

Only sequential implementation is discussed here. The present data structure is intrinsically sequential. In 
parallel implementations, a central time-histogram can distribute the unit next-time and state update calculation 
on several processors, with the drawback that the slower calculation limits the performance. Another idea is to 
consider several time-histograms on different processors and synchronize them. See (Mouraud et al., 2006) and 
(Morrison et al., 2007) for developments of these ideas. 

' Source code available at http: //enas . gf orge . inria . f r. 
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The fact that we use such a tiny simulation kernel has several practical advantages, e.g. to use spiking network 
mechanisms in embedded systems, etc.. However, it is clear that this is not "yet another" simulator because a com- 
plete simulator requires much more that an event queue (Brette et al., 2007). On the contrary, the implementation 
has been designed to be used as a plug-in in existing simulators, mainly MVASpike (Rochel & Martinez, 2003). 

2.2 From next event-time to lower-bound estimation 

Let us now consider the following modification in the event-based simulation paradigm. Each unit provides: 

-1'- the calculation of either the next event-time, 
or the calculation of a lower-bound of the next event-time, 

-2- the update of the neural unit when an internal or external event occurs, 
with the indication whether the previously provided next event-time or lower-bound is still valid. 

At the network level the mechanism's loop is now: 

-a' - retrieve either the next event-time and proceed to -b'- or a lower-bound and proceed to -c- 
-b'- require the update of the state of this unit, inform efferent units that this unit has emitted an 
event, and update the related event-times only if this event-time is lower than its previous estimation, 
-c- store the event-time lower-bound in order to re-ask the unit at that time. 

A simple way to interpret this modification is to consider that a unit can generate "silent events" which write: 
"Ask me again at that time, I will better know". 

As soon as each unit is able to provide the next event-time after a finite number of lower-bound estimations, the 
previous process is valid. 

This new paradigm is fully compatible with the original, in the sense that units simulated by the original 
mechanism are obviously simulated here, they simply never return lower-bounds. 

It appears that the implementation of this "variant" in the simulation kernel is no more than a few additional 
Une of codes. However, the specifications of an event-unit deeply change, since the underlying calculations can be 
totally different. 

We refer to this modified paradigm as the lazy event-based simulation mechanism. 
The reason of this change of paradigm is twofold: 

• Event-based and clock-based calculations can be easily mixed using the same mechanism. 

Units that must be simulated with a clock-based mechanism simply return the next clock-tick as lower- 
bound, unless they are ready to fire an event. However in this case, each unit can choose its own clock, 
requires low-rate update when its state is stable or require higher-rate update when in transient stage, etc.. 
Units with different strategies can be mixed, etc.. 

For instance, in (Morrison et al., 2005) units corresponding to synapses are calculated in event-based mode, 
while units corresponding to the neuron body are calculated in clock-based mode, minimizing the overall 
computation load. They however use a more complicated specific mechanism and introduce approximations 
on the next spike- time calculations. 

At the applicative level, this changes the point of view with respect to the choice of event-based/clock- 
based simulations. Now, an event-based mechanism can always simulate clock-based mechanisms, using 
this useful trick. 



6 



• Computation time can be saved by postponing some calculations. 

Event-based calculation is considered as costless than clock-based calculation because the neuron state is not 
recalculated at each time-step, only when a new event is input. However, as pointed out by several authors, 
when a large amount of events arrive to a unit, the next event-time is recalculated a large amount of time 
which can be much higher than a reasonable clock rate, inducing a fall of performances. 

Here this drawback can be limited. When a unit receives an event, it does not need to recalculate the next 
event-time, as soon as it is known as lower that the last provided event-time bound. This means that if the 
input event is "inhibitory" (i.e. tends to increase the next event-time) or if the unit is not "hyper-polarized" 
(i.e. not close to the firing threshold, which is not trivial to determine) the calculation can be avoided, while 
the opportunity to update the unit state again later is to be required instead. 

Remarks 

Mixing event-based and clock-based calculations that way is reasonable, only because the event-time queue 
retrieve/update operations have a very low cost. Otherwise, clock-ticks would have generated prohibitory time 
overloads. 

Changing the event-based paradigm is not a simple trick and may require to reconsider the simulation of neural 
units. This is addressed in the sequel for two important classes of biologically plausible neural units, at the edge of 
the state of the art and of common use: Synaptic conductance based models (Destexhe, 1997) and spike response 
models (Gerstner & Kistler, 2002). 



3 Event-based simulation of adaptive non-linear gIF models 

Let us consider a normalized and reduced "punctual conductance based generalized integrate and fire" (gIF) neural 
unit model (Destexhe, 1997) as reviewed in (Rudolph & Destexhe, 2006). 

The model is normalized in the sense that variables have been scaled and redundant constants reduced. This is 
a standard usual one-to-one transformation, discussed in the next subsection. 

The model is reduced in the sense that both adaptive currents and non-linear ionic currents are no more explic- 
itly depending on the potential membrane, but on time and previous spikes only. This is a non-standard approxi- 
mation and a choice of representation carefully discussed in appendix 3.1. 

Let V be the normalized membrane potential and (it = {• • • • • • } the list of all spike times < t. Here is 

the n-th spike-time of the neuron of index i. The dynamic of the integrate regime writes: 

— +g{t,uJt)v = i(t,LOt), (1) 

while the fire regime (spike emission) writes v(t) = 1 v{t'^) = with a firing threshold at 1 and a reset 
potential at 0, for a normalized potential. 
Equation (1) expands: 

^ + 1 _ E^] + ^ ^ (i _ t-) [V - E,\ = UQ,), (2) 

i n 

where tl and Ei are the membrane leak time-constant and reverse potential, while rj () and Ej the spike responses 
and reverse potentials for excitatory/inhibitory synapses and gap-junctions as made explicit in appendix A. Here, 
imO is the reduced membrane current, including simpUfied adaptive and non-Unear ionic current. 
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3.1 Reduction of internal currents 

Let us now discuss choices of modeling for im = /'"'"p* + 

Adaptive current 

In the Fitzhugh-Nagumo reduction of the original Hodgkin-Huxley model (Hodgkin & Huxley, 1952) the average 
kinematics of the membrane channels is simulated by a unique adaptive current. Its dynamics is thus defined, 
between two spikes, by a second equation of the form: 

JTadp 

^ = ~ ~ ^ + A. 5 (y - ^threshold) > 

with a slow time-constant adaptation ~ 144ms, a sub- threshold equivalent conductance (?^ ~ AnS and a 
level ~ O.OOSnA of spike-triggered adaptation current. It has been shown (Izhikevich, 2003) that when a 
model with a quadratic non-linear response is increased by this adaptation current, it can be tuned to reproduce 
qualitatively all major classes of neuronal in-vitro electro-physiologically defined regimes. 
Let us write: 

I-''PiV,t) = 7«*(to) + ^ /,* e^{V{s) - EL)ds + #{to,t) 

~ r'^Pito) + (l - e^^) {V - El) + A^ #{to,t) 

, ... spike— time dependent 

slow variation 

where #{to, t) is the number of spikes in the [to, t] interval while V is the average value between t and to, the 
previous spike-time. 

Since the time-constant adaptation is slow, and since the past dependency in the exact membrane potential value 
is removed when resetting, the slow- variation term is almost constant. This adaptive term is mainly governed by 
the spike-triggered adaptation current, the other part of the adaptive current being a standard leak. This is also 
verified, by considering the hnear part of the differential system of two equations in V and I°-'^p, for an average 
value of the conductance (5+ ~ 0.3 •• • 1.5n5 and G~ ~ 0.6 • • • 2.5715. It appears that the solutions are defined 
by two decreasing exponential profiles with ti ~ 16ms << T2 — 115ms time-constants, the former being very 
close to the membrane leak time-constant and the latter inducing very slow variations. 

In other words current adaptation is, in this context, mainly due to spike occurrences and the adaptive current 
is no more directiy a function of the membrane potential but function of the spikes only. 



Non-linear ionic currents 

Let us now consider the non-linear active (mainly Sodium and Potassium) currents responsible for the spike gen- 
eration. In models designed to simphfy the complex structure of Hodgkin-Huxley equations, the sub-threshold 
membrane potential is defined by a supra-linear kinematics, taken as e.g. quadratic or exponential, the latter form 
closer to observed biological data (Brette & Gerstner, 2005). It writes, for example: 



with Ea — — 40my the threshold membrane state at which the slope of the 1-V curve vanishes, while 5a = 2mV 
is the slope factor which determines the sharpness of the threshold. There is no need to define a precise threshold 
in this case, since the neuron fires when the potential diverges to infinity. 



>(y) = £i^e^>Owith 
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A recent contribution (Touboul, 2008) re-analyzes such non-linear currents, proposes an original form of the 
ionic current, with an important sub- threshold characteristic not present in previous models (Izhikevich, 2003; 
Brette & Gerstner, 2005) and show that one obtains the correct dynamics, provided that the profile is mainly 
non-negative and strictly convex. This is not necessarily a quadratic or exponential function. 

Making profit of this general remark, we propose to use a profile of the form of (3), but simply freeze the value 
of V to the the previous value obtained at the last spike time occurrence. This allows us to consider a supra-Unear 
profile which depends only on the previous spike times^. This approximation may slightly underestimate the ionic 
current before a spike, since V is increasing with time. However, when many spikes are input, as it is the case for 
cortical neurons, errors are minimized since the ionic current update is made at high rate. 

At a phenomenological level (Izhikevich, 2003) the real goal of this non-linear current in synergy with adaptive 
currents is to provide several firing regimes. We are going to verify experimentally that even coarser approxima- 
tions allow to attain this goal. 



3.2 Derivation of a spike-time lower-bound. 

Knowing the membrane potential at time to and the list of spike times arrival, one can obtain the membrane 
potential at time t, from (1): 

v{t) = u{to,t,u)to)vito)+ i/(s,i,a;tji(s,a;tjrfs (4) 

Jta 



^In fact a more rigorous result can be derived, altliougli at tlie implementation level, the simple heuristic proposed here seems sufficient. Let 
us write i(V, t, uJt) = i'{t. Cut) + /('°") (V, t, tit) thus separate the from all other currents written i' [t, tit). Let us consider the last 

spike time to of this neuron and let us write V the solution of the linear differential equation "withouf ' the ionic current ; 

C^+g{t,u,t)V = i'{t,u,t) 

with V{to) = Vreset, as obtained above. Define now V = V — V, with V{to) = 0, V being the solution of the previous equation (without 
the ionic current). This yields: 

C ^ + g(L Cjt)V = 7('°") (V + V{t, u:t), t, cDt) 
as easily obtained by superposition of the hnear parts of the equation. 

Let h(t, lit) be any regular function and f{V) any bijective regular function with f{V) ^ 0. These two functions allow us to model a whole 
family of ionic currents: 

liion) + V^t,U,t),t,U,t)= 9{t, U,t)V+ ^'jf^ ■ 

The choice of h and / is simply related to specific properties: The reader can easily verify, by a simple integration, that it allows to obtain a 
closed form: 

V{t, lit) = F-^ {^Jl^ h{s, Cjt) ds^ with F' = f and F(0) = 0. 

so that V is now a function of wt,t with V{to,Cdt) = 0, and so is /('°") {V{t, oJt) + V{t, a)t), t, a;*), removing the direct dependence on V. 
In other words, it now depends only on t and on the spike times (thus on uJt) and not anymore on the membrane potential explicitiy. Clearly, 
this only applies to neurons which have fired at least once during the period of observation. Otherwise, we assume that its initial condition was 

also Vreset- 

We can, .e.g, choose: 

Eait,C;jt) = V(t,u>t)-Saln[^^^if^) 
for any g > which allows to control the threshold for different conductance. 
Here h = g and f(v) = (ke^a — v)^^ for some k. 

In this case the threshold is no more fixed, but adaptive with respect to g{t): the higher the conductance, the higher the threshold (via the V). 
This is coherent with what has been observed experimentally (Azouz & Gray, 2000; Wilson, Weyrick, N. E. HaUworth, & Bevan, 2004), since 
the higher the conductance, the higher the frame rate increases with the spiking threshold. 
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with: ^ 

log{p{to,ti,u)to)) = - 9{s,(^to)ds (5) 

J to 

Furthermore, 

9{t,^to) = ^ + E,E„r.(i-i,") > 

iit,u;to) = ^EL + EjEnri{t-t]) Ej+irr^iQto) > 

since the leak time-constant, the conductance spike responses are positive, the reverse potential are positive (i.e. 
they are larger than or equal to the reset potential) and the membrane current is chosen positive. 

The spike-response profile schematized in Fig. 1 and a few elementary algebra yields to the following bounds: 

t e [to,ti] r^{to,h) < r{t) < ry{to,ti)+tr'y{toM) (7) 



writing r' the time derivative of r, with 

r/\{to,ti) = min(r(to),?'(ti)) andrv(to,ii) = max(r(to), r(ti)) 
with a similar definition for r(^(to, ^i)- 

Here we thus consider a constant lower-bound Ta and a Unear or constant upper-bound ry -\-try. The related 
two parameters are obtained considering in sequence the following cases: 





to e 


ti 6 


rv(to,ti) 


r'v{to,ti) 


(i) 


[*l - r(ii)/r'(ti),ta] 




r(ti)-fir'(fi) 


r'(ti) 


(ii) 


'/„./,,] 




r{h,) - lu r'(lu) 


/•'(/„) 


(111) 


l/...— ^1 


](„. +;x.^ 


(/■(/„) Li - r(ii) (o)/(^l - k,) 


('■I'lJ - '■(lu))/(l.i - lu) 


(iv) 


J - 'OC., <il 




Ktl) 





(V) 


J — oo, tbl 


[tb, +ooJ 


r{tb) 





(vi) 


[tb.flj 


[to, +oo\ 


r(to) 






In words, conditions (i) and (ii) correspond to the fact that the the convex profile is below its tangent (schema- 
tized by d in Fig. 1), condition (iii) that the profile is concave (schematized by d" in Fig. 1). In other cases, it can 
be observed that 1st order (i.e. linear) bounds are not possible. We thus use constant bounds (schematized by d' in 
Fig. 1). Conditions (iv) and (vi) correspond to the fact the profile is monotonic, while condition (v) corresponds to 
the fact the profile is convex. Conditions (iv), (v) and (vi) correspond to all possible cases. Similarly, the constant 
lower bound corresponds to the fact the profile is either monotonic or convex. 

From these bounds we derive: 



9{t,^to) > 

def 



■l"A(to,tl) 

=^ TeEl + Ej E„ ''i V (to , * 1 ) Ej + im (Wto ) 
+ i [E,En^^v(^0,tl)i?, 

while i'yiU); ti) > as the positive sum of r'^y values is always positive in our case. 
Combining with (4), since values are positive, yields: 



(8) 



v{t) < vy{t) =^ {v{to) - Vo) e 



(9) 
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Figure 1: The spike response profile. It has a flat response during the absolute delay interval [0,ta], an increasing convex 
profile until reaching its maximum at tb, followed by a decreasing convex and then concave profile, with an inflexion point at 
tc. After td the response is negligible. See text for details about d, d' and d". 



writing: 



def 

V, = T/\{to,ti) {iw{to,ti) - T/y{to,ti)i'^{to,ti)) 

def , . , 

Vo = V,+ I, to 

def 



(10) 



Finally we can solve the equation for = ^ (1) and obtain: 

' to 



io + rAlog(^^) 



vitO) > 1 
i'v(io,ii) > 
i'v{to,ti) = 0,i;. > 1 
otherwise 



(11) 



Here y = L{x) is the Lambert function defined as the solution, analytic at 0, of y = x and is easily tabulated. 

The derivations details are omitted since they have been easily obtained using a symbolic calculator. 

In the case where i[y{to, ti) = thus i, = 0, V\/{t) corresponds to a simple leaky integrate and fire neuron 
(LIF) dynamics and the method thus consists of upper bounding the gIF dynamics by a LIF in order to estimate 
a spiking time lower-bound. This occurs when constant upper-bounds is used for the currents. Otherwise (9) 
and (10) corresponds to more general dynamics. 



3.3 Event-based iterative solution 

Let us apply the previous derivation to the calculation of the next spike-time lower-bound for a glF model, up to a 
precision Cf. One sample run is shown Fig. 2 

Given a set of spike-times and an interval [ta, ti], from (8) Ta, ?'v and I'y are calculated in about 10^ M C 
operations, for an average of C connections per units, and with an average number M < N of firing units. This 
is the costly part of the calculation^, and is equivalent to a single clock-based integration step. Spike response 

'it appears, that since each synapse response corresponds to a linear differential equation which could be analytically integrated, the global 
synaptic response fj{t) = rj{t — t^) can be put in closed form, and then bounded by constant values, thus reducing the computation 
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Figure 2: An example of gIF normalized membrane potential. The left trace corresponds to 50ms of simulation, 
the neuron being artificially connected to a periodic excitatory/inhibitory input neuron pair at 50/33 Hz of high 
synaptic weight, in order to make exphcit the double exponential profiles. The right trace corresponds to 200ms 
of simulation, with a higher inhibition. The weights have been chosen, in order to make the neuron adaptation 
expUcit: the firing frequency decreases until obtaining a sub-threshold membrane potential. 



profiles rj() and profile derivatives r^O for excitatory/inhibitory synapses and gap junctions are tabulated with a 
et step. Then, from (10) and (11) we obtain t'^{to, ti). 

The potential v{to) is calculated using any well-estabUshed method, as detailed in, e.g., (Rotter & Diesmatm, 
1999), and not reviewed here. 

The following algorithm guarantees the estimation of a next spike-time lower-bound after to. Let us consider 
an initial interval of estimation d, say d 10ms: 



-a- The lower-bound ty = (^o ■ ^^o + is calculated. 

the lower-bound value to + dis returned 
the estimation interval is doubled d •*— 2 

the lower-bound value ty is returned 
the estimation interval is reduced d ^ 
-d- If to < ty < to + d, d < et, the next spike- time is returned. 



-h-lfto + d< ty 
-c- If to + et < *v < ^0 + 



l/^/2d 



Step -b- corresponds to the case where the neuron is not firing in the estimation interval. Since v{t) is bounded 
by Vy{t) and the latter is reaching the threshold only outside [to, to + d\,to + d\s,a. time lower-bound. In addition, 
a heuristic is introduced, in order to increase the estimation interval in order to save computation steps. 

Step -c- corresponds to a strict lower-bound computation, with a relative value higher than the precision . 

Step -d- assumes that the lower-bound estimation converges towards the true spike-time estimation when ti 

to. 

This additional convergence property is easy to derive. Since: 

limti^to rjA{to,ti) = limt^^to rjv{to,ti) = rj{to),limti^to r'jy{to,ti) = rj{to), 

then: 

limti^to l/T/^{to,ti) = fif(to),limti^to iw{to,ti) = i(to), limtj^to «'v(*o,*i) = i'{to), 



complexity from 0{MC) to 0(C) as detailed in (Rotter & Diesmami, 1999). This weU-known issue is not re-addressed here, simply to avoid 
making derivations too heavy. 
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yielding: 

limt,^t„ vv{t) = v^{t) =' {v{to) - Vo) e-9(*°) +v,+i,t 

writing: 

i. l/g{to)i'{to) 

V, l/9ito){i{to)-l/g{to)i'{to)) 

dcf _ , ^ , 
Vo = V, +t, to 

We thus obtain a limit expression V'^ of t^v when ti ^ to. From this limit expression we easily derive: 

v{t) - v^{t) = -1/2 g' {to) v{to) f + 0{t% 
and finally obtain a quadratic convergence, with an error closed-form estimation. 

The methods thus corresponds to a semi-interval estimation methods of the next spike-time, the precision 
being adjustable at will. 

The interval of estimation d is adjusted by a very simple heuristic, which is of standard use in non-linear 
numerical calculation adjustment. 

The unit calculation corresponds to one step of the iterative estimation, the estimation loop being embedded 
in the simulator interactions. This is an important property as far as real-time computation is concerned, since a 
minimal amount of calculation is produced to provide, as soon as possible, as suboptimal answer. 

This "lazy" evaluation method is to be completed by other heuristics: 

- if the input spike is inhibitory thus only delaying the next spike-time, re-calculation can be avoided; 

- if all excitatory contributions g\/ — M (tb) are below the spiking threshold, spike-time is obviously infinity; 

- after a spike the refractory period allows us to postpone all calculation (although synaptic conductances still 
integrate incoming spikes). 

In any case, comparing to other glF models event-based simulation methods (Brette, 2006, 2007; Rudolph & 
Destexhe, 2006), this alternative method allows one to control the spike-time precision, does not constraint the 
synaptic response profile and seems of rather low computational cost, due to the "lazy" evaluation mechanisms. 

Numerical convergence of the lower-bound iteration 

Considering biologically plausible parameters as reviewed here and in appendix A, we have experimented carefully 
the numerical convergence of this lower-bound iterative estimation, considering a glF neuron with adaptive and 
non-linear internal currents, implemented as proposed here, and providing membrane potential, e.g., as shown in 
Fig. 2. 

Let us report our numerical experimentation. 

We have always observed the convergence of the method (also extensively experimented at the network level in 
a next section), with a convergence in about 2—20 iterations (mean ~ 11, standard-deviation ~ 5), the lower-bound 
iteration generating steps of about 0.01 — lOm.s (mean c± 3ms, standard-deviation ~ 4ms) from one lower-bound 
to another, with three distinct qualitative behaviors: 

-a- Sub-threshold maximal potential: the previous calculation estimates a maximal membrane potential below the 
threshold and calculation stops, the neuron being quiet; in this mode the event-based strategy is optimal and a large 
number of calculations are avoided with respect to clock-based paradigms. 

-b- Sub-threshold lower-bound estimation: the maximal membrane potential is still estimated over the threshold, 
with a next-spike time lower bound. In this mode, we observed an exponential increase of the next-spike time 
lower bound and in 2 — 5 iterations the maximal membrane potential estimation is estimated under the threshold, 
switching to mode -a-; in this mode the estimation interval heuristic introduced previously is essential and the 
next-spike time lower bound estimation allows the calculation to quickly detect if the neuron is quiet, 
-c- Iterative next-spike time estimation: if a spike is pending, the previous calculation estimates in about 10 — 20 
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iterations the next-spike time up to a tunable precision (corresponding to the dt of the simulation mechanism). The 
present mechanism acts as an iterative estimation of the next spike time, as expected. 

4 About event-based simulation of SRM models 

Among spiking-neuron models, the Gerstner and Kistler spike response model (SRM) (Gerstner & Kistler, 2002) 
of a biological neuron defines the state of a neuron via a single variable: 

where Ui is the normalized membrane potential, j() is the continuous input current for an input resistance ri. The 
neuron fires when Ui{t) > 6i, for a given threshold, Ui describes the neuronal response to its own spike (neuronal 
refractoriness), t* is the last spiking time of the ith neuron, e,, is the synaptic response to pre-synaptic spikes at 
time fj post-synaptic potential (see Fig. 3), Wij is the connection strength (excitatory if Wij > or inhibitory if 
Wij < 0) and 5ij is the connection delay (including axonal delay). Here we consider only the last spiking time t* 
for the sake of simplicity, whereas the present implementation is easily generalizable to the case where several are 
taken into account. 




Figure 3: Potential profiles used in equation (12). The original exponential profiles derived by the authors of the 
model are shown as thin curves and the piece-wise linear approximation as thick Unes. Any other piece- wise linear 
profiles could be considered, including piece- wise finer linear approximation of exponential profiles. 

Let us call here ISRM such piece- wise linear approximations of SRM models. 

This model is very useful both at the theoretical and simulation levels. At a computational level, it has been 
used (see (Maass & Bishop, 2003) for a review) to show that any feed-forward or recurrent (multi-layer) analog 
neuronal network (a-la Hopfield, e.g., McCuUoch-Pitts) can be simulated arbitrarily closely by an insignificantly 
larger network of spiking neurons, even in the presence of noise, while the reverse is not true (Maass, 1997; Maass 
& Natschlager, 1997). In this case, inputs and outputs are encoded by temporal delays of spikes. These results 
highly motivate the use of spiking neural networks. 

This ISRM model has also been used elsewhere (see e.g. (Schrauwen, 2007) for a review), including high-level 
specifications of neural network processing related to variational approaches (Vieville, Chemla, & Kornprobst, 
2007), using spiking networks (Vieville & Rochel, 2006). The authors used again a ISRM to implement their 
non-linear computations. 

Let us make expUcit here the fact a ISRM can be simulated on an event-based simulator for two simple reasons: 

- the membrane potential is a piece-wise linear function as the sum of piece-wise linear functions (as soon as the 
optional input current is also piece-wise constant or linear), 

- the next spike-time calculation is obvious to calculate on a piece-wise linear potential profile, scanning the linear 
segments and detecting the 1st intersection with u = 1 if any. The related piece- wise linear curve data- structure 
has been implemented' and support three main operations: 
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- At each spike occurrence, add linear pieces to the curve corresponding to refractoriness or synaptic response. 

- Reset the curve, after a spike occurrence. 

- Solve the next spike-time calculation. 
This is illustrated in Fig. 4. 

SRH trace 1 SRH trace 2 

1 I 1 1 1 1 1 1 I 1 1 1 1 1 




Figure 4: An example of ISRM normalized membrane potential. The traces corresponds to 100ms of simulation. 
The leftward trace uses the fastest possible piece-wise linear approximation of the SRM model profiles. The 
rightward trace is simulated with a lower excitation inhibition and using a thinner piece- wise linear approximation. 
The neuron is defined by biologically plausible parameters as reviewed in appendix A . It is coimected to a pair of 
periodic excitatory and inhibitory input neurons, with different time constants, as in Fig. 2. 

This is to be compared with other simulations (e.g. (Mattia & Giudice, 2000; Ruf, 1998)) where stronger 
simplifications of the SRM models have been introduced to obtain a similar efficiency, whereas other authors 
propose heavy numerical resolutions at each step. When switching from piece-wise Unear profiles to the original 
exponential profiles (Gerstner & Kistler, 2002), the equation to solve is now of the formal form: 

without any closed-form solution as soon as n > 1. One elegant solution (Brette, 2007) is to approximate the 
time-constant by rational numbers and reduce this problem to a polynomial root finding problem. Another solution 
is to upper-bound the exponential profiles by piece-wise linear profiles in order to obtain a lower-bound estimation 
of the next spike-time and refine in the same way as what has been proposed in the previous section. Since the 
mechanism is identical, we are not going to further develop. 

In any case, this very powerful phenomenological model of biological dynamics can be simulated with several 
event-based methods, including at a fine degree of precision, using more complex piece-linear profiles. 

5 Experimental results 

5.1 Kernel performances and features 

In order to estimate the kernel sampling capability we have used, as a first test, a random spiking network with 
parameter less connections, the spiking being purely random thus not dependent on any input. 

In term of performances, on a standard portable computer (Pentium M 750 1.86 GHz, 512Mo of memory) we 
process about 10^"*^ spike-time updates / second, given the network size and coimectivity. Performances reported 
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in Fig. 5 confirm that the algorithmic complexity only marginally depends on the network size, while it is mainly 
function of the number of synapses (although both quantities are indeed linked). We also notice the expected tiny 
overhead when iterating on empty boxes in the histogram, mainly visible when the number of spikes is small. 
This overhead is constant for a given simulation time. The lack of proportionality in performances is due to the 
introduction of some optimization in the evaluation of spike-times, which are not updated if unchanged. 

We have also observed that the spike-time structure upper and lower bounds D and dt have only a marginal 
influence on the performances, as expected. 

Moreover, the numbers in Fig. 5 allows to derive an important number: the overhead for an event-based imple- 
mentation of a clock-based mechanism. Since we can process about 2 10® updates/second while we have measured 
independently that a minimal clock-based mechanism process about 5 10^ updates/second, we see that the cost the 
overhead is of about 0.5/is / update. This number is in coherence with the number of operations to realize at time 
modification in the underlying data structure. 

kernel perf ornances 

7 I 1 1 1 1 1 1 




g I 1 1 , , 1 1 

3 3,3 A 4.3 3 3.5 6 

size <L0{<(1G>> 

Figure 5: Simulation performances, for 10^ spike firings and about 0.1s of simulation time. The CPU time is 
about 2s for 2^^ neurons and 2^^ synapses and does not depend on D or dt values, as expected. Spike-time 
insertion/deletion are counted as elementary updates. The network size in abscissa varies from 10'^ to about 10® 
and the number of connections from from to 10^, corresponding to curve end-points. Curves are shown for a 
connection probability of P = (black, upper curve), P — 10^"^ (brown), P = 10^^ (red), P = 10^^ (orange). 
The performance is mainly function of the number of synapses, and marginally of the number of neurons. 

It is important to clarify these apparently "huge" performances. The reason is that the event-based simulation 
kernel is minimal. As detailed in Table 1, the implementation make a simple but extensive use of the best mecha- 



16 



nisms of object oriented implementations. The network mechanism (i.e., the kernel) corresponds to about 10Kb of 
C++ source code, using a 0{D/dt + N) buffer size and about 0(1 + C + e/dt) ~ 10 — 50 operations/spike, for 
a size N network with C connections in average, while e <C 1. This e corresponds to the overhead when iterating 
on empty boxes in the histogram. We can use such a simple spike-time data structure because of the temporal 
constraints taken into account in our specifications. 

As a consequence, not all spiking mechanisms are going to be simulated with this kernel: units with event-time 
intervals or input/output event delay below dt are going to generate a fatal error; units with inter-event intervals 
higher than D are also going to defeat this mechanism (unless an extension of the present mechanism, discussed 
previously, is not implemented). Note that despite these limitations the event-time accuracy itself is the not dt but 
the floating point machine precision. 

template <class C> class Unit { 

// Gets the next alter time: event time or a lower-bound of the next event time, 
inline virtual double getNext (double present-time) ; 

// Called to update the unit state, when the next spiking time or a lower-bound occurs. 

// Returns true if an event occurs, false it was a lower-bound 
inline virtual bool next (double present-time); 

// Called when an input unit fires an event. 

// Returns true if the next alert time is to be updated, false otherwise. 

inline virtual bool add(int neuron-index, C& connection-parameter, double present-time); 

}; 

Table 1: Specification of an event-based unit (pseudo-code). Each unit (neuron or group of neurons) specifies is 
next "alert" time and informs the network about event-occurrence. Lazy evaluation is implemented, at this level, 
via the fact that alert time is optional updated when receiving an event. The connection is templated in order for 
the kernel to be optimally recompiled for each kind of connection, while unit's mechanisms are inlined, allowing 
the compiler to eUminate code interface. The connection parameters is passed by reference, in order adaptation 
mechanisms to be implemented. See text for further details. 



Clock-based sampling in event-based environment. A step further, we have implemented a discretized version 
of a glF network, called BMS, as detailed in (Cessac, 2008; Cessac & Vieville, 2008) (equations not reported here). 
The interest of this test is the fact that we can compare spike by spike an event-based and clock-based simulation 
since the latter is well-defined, thus without any approximation with respect to the former (see (Cessac & Vieville, 
2008) for details). 

We have run simulation with fully connected networks of size, e.g. N = 100 — 1000 over observation periods 
of T = 1000 — 10000 clocks, with the same random initial conditions and the same randomly chosen weights, as 
show in Fig.6. We have observed: 

-1- the same raster (i.e., with a Victor-Purpura distance of (Victor, 2005)); this exactitude is not surprising 
despite the fact that floating point errors accumulate^: we are performing the same floating point calculations in 
both cases, since the event-based implementation is exact, thus . . with the same errors; 

-2- the overhead of the event-based implementation of the clock-based sampling is negligible (we obtain a 
number < Q.ljis/ step), as expected. Again this surprisingly slow number is simply due to the minimal implemen- 
tation, based on global time constraints, and the extensive use of the C/C++ optimization mechanisms. 

'*Note that even if time is discretized, for BMS networks, the dynamics is based on floating point calculations, thus floating point errors 
accumulate. However as soon as spike is fired, the potential is reset and previous errors are canceled. This explains why time-discretized 
simulations of IF networks are numerically rather stable. 
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Figure 6: An example of BMS neural network simulation used to evaluate the clock-based sampling in the event- 
based kernel, for TV = 1000 and T = 10000 thus 10^ events. The first 100 neurons activity is shown at the end of 
the simulation. The figure simply shows the strong network activity with the chosen parameters, and . 



Kernel usage. A large set of research groups in the field have identified what are the required features for such 
simulation tools (Brette et al., 2007). Although the present implementation is not a simulator, but a simple simula- 
tor plug-in, we have made the exercise to list to which extends what is proposed here fits with existing requirements, 
as detailed in Table 3 . We emphasize the fact the required programming is very light, for instance a "clock" neuron 
(allowing to mix clock-based and event-based mechanisms) writes: 

class ClockUnit : public Unit<bool> { 

ClockUnit (double DT) : DT(DT), t (DT) double DT, t; 

inline virtual double getNext (double present-time) return t; ; 

inline virtual bool next (double present-time) t += DT; return true; ; 

inline virtual bool add(int neuron-index, bool connection-parameter, double present-time) return 

false; ; 

}; 

providing DT is the sampling period. It is not easy to make things simple, and several possible choices of imple- 
mentations have been investigated before proposing the interface proposed in Table 1 . 

See (Rochel, 2004) for a further description of how event-based spiking neuron mechanisms can be imple- 
mented within such framework. Although presented here at a very pragmatical level, note that these mechanisms 
are based on the modular or hierarchical modeling strategy borrowed from the DEVS formalism (see, e.g., (Rochel, 



18 



2004)). 



5.2 Experimenting reduced adaptive and ionic currents 

In order to experiment about our proposal to reduce ionic and adaptive currents to a function depending only on 
spike time, we consider a very simple model whose evolution equation at time t for the membrane potential v is: 

If (t = 0) 

V = 0; u = 0; t.O = 0; 

else if (v > 1) 

V = 0; u = u + -k; t.O = t; (13) 

else 

V = -g (v - E) - u + i; 
if (t > t.O + d) V = 

where u is the adaptive current (entirely defined by equation (13)), t_0 the last spiking time, d the non-linear 
current delay. The differential equation is simulated using an Euler interpolation as in (Izhikevich, 2003; Touboul, 
2008) to compare our result to what has been obtained by the other authors. The obvious event-based simulation 
of this model has been also implemented'. The input current 1 is either a step or a ramp as detailed in Fig. 7 and 
Table 2. 

Four parameters, the constant leak conductance g, the reversal potential E, the adaptation current k step and 
the (eventually infinite) non-linear current delay d allows to fix the firing regime. These parameters are to be 
recalculated after the occurrence of each internal or external spike. In the present context, it was sufficient to use 
constant value except for one regime, as made expUcit in Fig. 2. We use the two-stages current whose action is to 
reset the membrane current after a certain delay. We made this choice because it was the simplest and leads to a 
very fast implementation. 

Experimental results are given in Fig. 7 for the parameters listed in Fig. 2. These results correspond to almost 
all well-defined regimes proposed in (Izhikevich, 2003). The parameter adjustment is very easy, we in fact use 
parameters given in (Izhikevich, 2003) with some tiny adjustments. It is an interesting numerical result: The 
different regimes are generated by parameters values closed to those chosen for the quadratic model, the dynamic 
phase diagram being likely similar. See (Touboul, 2008) for a theoretic discussion. 

This places a new point on the performance/efficiency plane proposed by (Izhikevich, 2004) at a very challeng- 
ing place, and see that we can easily simulate different neuronal regimes with event-based simulations. 

However, it is clear that such model does not simulates properly the neuron membrane potential as it is the 
case for the exponential model (Brette & Gerstner, 2005). It is usable if and only if spike emission is considered, 
whereas the membrane potential value is ignored. 
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Table 2: Examples of parameters used to generate the spiking modes shown in Fig. 7. The mixed mode is simulated 
by a variable adaptation step fc = {— 20, 20}. 
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Figure 7: Typical resuhs showing the versatihty of the reduced model for spiking, bursting and other modes, includ- 
ing and different current-frequency-responses (CFR). For each mode, the upper trace shows the action potentials, 
the lower trace the input current. These results include the excitatory mode of type I where the spike frequency 
can be as small as possible in a 1 — 10'^ Hz range and of type II where the spike frequency remains bounded. Tonic 
spiking is not shown, since obvious to obtain. 



5.3 Experimental benchmarks 

We have reproduced the benchmark 4 proposed in Appendix 2 of (Brette et al., 2007) which is dedicated to 
event-based simulation: it consists of 4000 IF neurons, which 80/20% of excitatory /inhibitory neurons, connected 
randomly using a connection probability of 1/32 ~ 3%. So called "voltage-jump" synaptic interactions are used: 
the membrane potential is abruptly increased/decreased by a value of 0.25/2.25mV^ for each excitatory/inhibitory 
event (thus using fixed randomly chosen weights). Here, we also introduce a synaptic delay of 2 /4ms respectively 
and an absolute refractory period of 1ms, both delays being corrupted by an additive random noise of lO^s of 
magnitude. We also have increased the network size and decreased the connection probability to study the related 
performances. In this network a synapse is simply defined by an index, weights are constant. See (Brette et al., 
2007) for further details. One result is proposed in Fig. 8 to qualitatively verify the related network dynamics. The 
fact we find small inter-spike intervals in this case is coherent with previous observed results. 

A step further, we also made profit of the new proposed approximation of gIF neuron models to run another 
test, inspired by another benchmark proposed in (Brette et al., 2007), after (Vogels & Abbott, 2005), considering 
current-based interactions (CUBA model) and/or conductance-based interactions (COBA model). In our context, 
current based interactions correspond to gap junctions, while conductance-based interactions correspond to synap- 
tic junctions. It was useless to reproduce the original benchmarks in (Brette et al., 2007) or (Vogels & Abbott, 
2005), but interesting to experiment if we can explore the network dynamics with the improved model proposed 
here, using the method proposed in 3.1, and the parameters reviewed in appendix A, thus beyond CUBA/COBA 
models. 
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Figure 8: Inter-spike interval histogram (left view) in linear coordinates measured after t > 0.1s, and correspond- 
ing raster plot (right-view) during Is of simulation, for the benchmark 4 proposed in Appendix 2 of (Brette etal, 
2007). 

One result is shown in Figs. 9. Results are coherent with what is discussed in details in (Vogels & Abbott, 
2005), and in particular close to what has been reviewed in (Vogels, Rajan, & Abbott, 2005). This is clearly a 
preliminary test and the influence, on the network dynamics, of thus alternate model is out of the scope of the 
present work, and a perspective for a further study. 

6 Discussion 

Taking global temporal constraints into account, it has been possible to better understand, at the simulation levels 
to which extends spiking mechanisms are bounded and simplified. At this simulation level, the challenge is to gen- 
erate spike-trains corresponding to what is observed in biology or to what is required for computational processing, 
without the necessity to precisely reproduce the internal neuron state. This is a very important simplification when 
the goal is to switch from the neural scale to the network one. 

The proposed mechanism is a complement of existing simulation tools (Brette et al., 2007) in the following 
quantitative and quahtative senses: 

Quantitative complementarity 

As a software module, it has been designed to be as fast as possible. 

The cost for this choice is that a programmatic interface is required, while in order to be available on any 
platform with the fastest performances, a C/C++ implementation is required (interfaces to other programming 
languages being available). The fact that it is also a "small kernel" allows us to target embedded applications: 
since computing with spikes is now a mature methodology, a tool to run such algorithms on various platforms 
(e.g., in robotics) or embedded systems (e.g., intelligent reactive devices) was required. 

This has been possible here, without any loss in precision, the underlying data-structures being strongly sim- 
plified, but with another drawback: the network dynamics is constrained since spiking units must verify temporal 
constraints. 

A step further, the use of models at the edge of the state of the art, such adaptive non-linear glF networks, or 
SRM networks is made possible in an event-based framework, thus with expected better performances. 

Regarding gIF networks, (Brette, 2007) has proposed a pure event-based method taking step-wise synapses 
with exponential decays into account. The same level of modeling has been proposed by (Morrison et al., 2005), 
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Figure 9: Inter-spike interval histogram for excitatory neurons (top-left) and inhibitory neurons (down-left) with 
the corresponding raster head (right). The abscissa is the decimal of the histogram log of the interval and the 
ordinate the inter-spike observed probability. 



mixed with clock-based mechanisms, while (Rudolph & Destexhe, 2006) have investigated how to take synaptic 
alpha profiles into account. The proposed methods are based on sophisticated analytical derivations with the risk of 
having a rather huge number of operations to perform at each step. As a complementary variant of these methods, 
we propose here to introduce another degree of freedom, using iterative lower-bound estimations of the next spike 
time. This heuristic applied to gIF neurons seems to converge quickly in practice. 

Regarding SRM networks, we have generalized the simple idea to use piece-wise linear profiles, approximating 
the original exponential profiles roughly (replacing the exponential curve by a line-segment) or at any level of 
precision (approximating the exponential curve by any number of line-segments). The precision/performance 
trade-off is thus adjustable at will. 

The reason to consider gIF and SRM neuron simulation here is that they correspond, up to our best knowledge, 
to the most interesting punctual neuron models actually used in biologically plausible neural network simulation, 
in the deterministic case. 

Qualitative complementarity 

Two key points allow us to performs new simulations with this tool: 

Event-based and clock-based mechanisms can be easily mixed here in an event-based simulation mechanism, 
whereas other implementations mix clock-based and event-based mechanisms in a clock-based simulator (e.g., in 
(Morrison et al., 2005)), or use spike-time interpolation mechanisms in order to better approximate event-based 
mechanisms in such clock-based environment. Using an event-based simulator to simulate a clock is obvious, but 
usually stupid, because the event-based mechanism usually generates an heavy overhead, thus making the clock- 
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based part of the simulation intractable. This is not the case here, since we use this minimal data-structure and 
have been able to see that the overhead is less than one micro-second on a standard laptop. It is thus appears a 
good design choice to use an event-based simulation mechanism to mixed both strategies. 

The second key point is that, we have proposed a way to consider adaptive non-linear gIF networks in an 
event-based framework. It is easy to get convinced that 2D integrate and fire neurons differential equations with 
non-Unear ionic currents (e.g. exponential, quartic or quadratic (Touboul, 2007)) do not have closed-form solutions 
(unless in very special cases). Therefore, the next spike time is not calculable, except numerically, and the exact 
event-based implementation is not possible. Alternatives strategies have been proposed such as simulations with 
constant voltage steps (Zheng, Tonnelier, & Martinez, 2008) allowing to implement quadratic ID (thus without 
adaptive currents) gIF networks in a modified event-base framework. In order to get rid of these limitations, one 
proposal developed here is to consider adaptive currents which depends only on the previous spiking time neuron 
state and non-linear ionic currents updated only at the each incoming spike occurrence. With these additional 
approximations, the event-based strategy can be used with such complex models. This is a complementary heuristic 
with respect to existing choices. 

We notice that the present study only considers deterministic models, while the simulation of stochastic net- 
works is also a key issue. Hopefully, event-based implementations of network of neurons with stochastic dynamics 
is a topic already investigated, both at the computer implementation level (Rochel, 2004) and modeling level 
(Reutimann, GuigUano, & Fusi, 2003). In the latter case, authors propose to reduce the multiple stochastic neu- 
ronal input activity to a dedicated stochastic input current, and investigate this choice of modeling in an event-based 
framework, making explicit very good performances. This method seems to be easily implementable in our present 
kernel, though this is out of the scope of the present work. 

Conclusion 

At a practical level, event-based simulation of spiking networks has been made available, using the simplest pos- 
sible programmatic interface, as detailed previously. The kernel usage has been carefully studied, following the 
analysis proposed in (Brette et al., 2007) and detailed in Table 3 

The present implementation thus offers a complementary alternative with respect to existing methods, and 
allows us to enrich the present spiking network simulation capabilities. 



A Appendix: About gIF model normalization 

Let us review how to derive an equation of the form of (2). We follow (Izhikevich, 2004; Brette & Gerstner, 2005; 
Rudolph & Destexhe, 2006) in this section. We consider here a voltage dynamics of the form: 

j jleak _|_ jsyn _|_ jgap jext _j_ jadp _|_ jion 

dt 

thus with leak, synaptic, gap-junction, external currents discussed in this section. 

Membrane voltage range and passive properties 

The membrane potential, outside spiking events, verifies: 

V{t) e [Vreseh ^threshold] 

with typically Vj-gset — — SOmV and a threshold value threshold — — 50mF± lOm V. When the threshold 
is reached an action potential of about 1-2 ms is issued and followed by refractory period of 2-4 ms (more precisely. 
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an absolute refractory period of 1-2 ms without any possibility of another spike occurrence followed by a relative 
refraction to other firing). Voltage peaks are at about 40mV and voltage undershoots about —90mV. The threshold 
is in fact not sharply defined. 

The reset value is typically fixed, whereas the firing threshold is inversely related to the rate of rise of the action 
potential upstroke (Azouz & Gray, 2000). Here it is taken as constant. This adaptive threshold diverging mecha- 
nism can be represented by a non-linear ionic current (Izhikevich, 2004; Brette & Gerstner, 2005), as discussed in 
section 3.1. 

From now on, we renormalize each voltage between [0, 1] writing: 

^ " Preset ^^^^ 



threshold ^eset 



The membrane leak time constant tl — 20ms is defined for a reversal potential El — —80 mV, as made 
explicit in (2). 

The membrane capacity C = S Cl ~ 300pF, where Cl — 1 ^iFcm~'^ and the membrane area S ~ 
38.013 /;im^, is integrated in the membrane time constant tl = Cl/Gl where Gl — 0.0452 mS'cm"-^ is the 
membrane passive conductance. 

From now on, we renormalize each current and conductance divided by the membrane capacity. Normalized 
conductance units are s^^ and normalized current units V/ s. 



Synaptic currents 

In conductance based model the occurrence of a post synaptic potential on a synapse results in a change of the 
conductance of the neuron. Consequently, it generates a current of the form: 

ry^{V,Cbt,t) = ^Gf{t,ijt) [V{t) - E+] + ^GT{t,Cjt) [V{t) - E_] , 

j 3 

for excitatory + and inhibitory — synapses, where conductances are positive and depend on previous spike-times 

In the absence of spike, the synaptic conductance vanishes (Koch, 1999), and spikes are considered having an 
additive effect: 

Gfit,u>t) = G^J2^^it^t]) 

n 

while the conductance time-course r^{t — f") is usually modeled as an "exponential", "alpha" (see Fig. 10) or 
two-states kinetic (see Fig. 11) profile, where H is the Heaviside function (related to causality). 

Note, that the conductances may depend on the whole past history of the network, via Cot. 

The "exponential" profile (r{t) = H{t) e~^) introduces a potentially spurious discontinuity at the origin. 
The "beta" profile is closer than the "alpha" profile to what is obtained from a bio-chemical model of a synapse. 
However, it is not clear whether the introduction of this additional degree of freedom is significant here. Anyway, 
any of these can be used for simulation with the proposed method, since their properties correspond to what is 
stated in Fig. 1. 

There are typically, in real neural networks, 10^ excitatory and about 2 10^ inhibitory synapses. The corre- 
sponding reversal potential are E^ ~ mV and E- ~ — 75 mV, usually related to AMPA and GABA receptors. 
In average: G^ ~ 0.66 nS',r+ ~ 2ms and (5j ~ Q.Q'^nS,T~ ~ 10 ms, for excitatory and inhibitory synapses 
respectively, thus about 570ms~^, 600ms~^ in normalized units, respectively. The coefficients (5^ give a measure 
of the synaptic strength (unit of charge) and vary from one synapse to another and are also subject to adaptation. 
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Figure 10: The "alpha" profile a{t) ^ H (t) ^ e ^ plotted here for r = 1. It is maximal at i = r with ^(t) = 1 /e, 

the slope at the origin is 1/r and its integral value /g^°° Oi{s) ds = t since (J a){t) = (t — <) e^^ + k. This 
profile is concave for t e]0, 2 t[ and convex for t G]2 r, +oo[, while a{2 r) ~ 2/e^ at the inflexion point. 




Figure 11: The two-states kinetic or "beta" profile /3{t) — Hit) -^^^ (e^^ ~ e^'^ is plotted with a normalized 
magnitude for the same r = 1 as the "alpha" profile but different k = 1.1, 1.5, 2, 5, 10 showing the effect of this 
additional degree of freedom, while \]m^^i[3{t) = a{t). The slope at the origin and the profile maximum can be 
adjusted independently with "beta" profiles. It is maximal att, = t ln{n) /{k — 1) the slope at the origin is l/r 
and its integral value t/k. The profile is concave for t ejO, 2 and convex for t e]2i,,+oo[. 

This framework affords straightforward extensions involving synaptic plasticity (e.g. STDP, adjusting the 
synaptic strength), not discussed here. 

Gap junctions 

It has been recently shown, that many local inter-neuronal connections in the cortex are realized though electri- 
cal gap junctions (Galarreta & Hestrin, 2001), this being predominant between cells of the same sub-population 
(Amitai et al., 2002). At a functional level they seem to have an important influence on the synchronicity between 
the neuron spikes (Lewis & Rinzel, 2003). Such junctions are also important in the retina (Wohrer & Kornprobst, 
2008). 

The electrotonic effect of both the sub-threshold and supra-threshold portion of the membrane potential Vj (t) 
of the pre-junction neuron of index j seems an important component of the electrical coupling. This writes (Wohrer 
& Kornprobst, 2008; Lewis & Rinzel, 2003): 

pap^y^ - G* [[V, it) - V{t)] + E, E„ rit - tp] 
where G* is the electrical coupling conductance, the term Vj{t) — V{t) accounts for the sub-threshold electrical 
influence while and E, parametrizes the spike supra-threshold voltage influence. 

Regarding the supra-threshold influence, a value E, ~ 80mV corresponds to the usual spike voltage magnitude 
of the spiking threshold, while r, ~ 1ms corresponds to the spike raise time. Here r() profile accounts for the 
action potential itself, slightly filtered by the biological media, while the gap junction intrinsic delay is of about 
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10/xs. These choices seem reasonable with respect to biological data (Galarreta & Hestrin, 2001; Lewis & Rinzel, 
2003) 

A step further, we propose to neglect the sub-threshold term for three main reasons. On one hand, obviously the 
supra-threshold mechanism has a higher magnitude than the sub-threshold mechanism, since related to action po- 
tentials. Furthermore, because of the media diffusion, slower mechanisms are smoothed by the diffusion, whereas 
faster mechanisms better propagate. On the other hand, this electrical influence remains local (quadratic decrease 
with the distance) and is predominant between cells of the same sub-population which are either synchronized or 
have a similar behavior, as a consequence \ Vj{t) — V{t)\ remains small for cells with non-negligible electrical 
coupUng. Furthermore, a careful analysis of such electrical coupling (Lewis & Rinzel, 2003) clearly shows that the 
sub-threshold part of the contribution has not an antagonist effect on the neuron synchrony, i.e., it could be omitted 
without changing qualitatively the gap junction function. 

As a conclusion, we are able to take gap junction into account with a minimal increase of complexity, since we 
obtain a form similar to synaptic currents, using very different parameters 

External currents 

Direct input (or external) current is often related to electro-physiological clamps. At another level of representation, 
the average activity of the neuron can be modeled by a constant or random input current. In both cases, the way the 
proposed simulation methods is proposed requires to assume such external current to be taken as constant between 
two spikes, i.e. having temporal variations small enough to be neglected. If not, it is easy to associate to the 
external current an event unit which fires at each new current value, in order the neuron to take this new value into 
account. 
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Features 


f~'lnpk'-V>a';prl 


■ pan It Qimiilfitp ploplf-fia^pH QtratptripQ 


yes 




■ iTi t n 1 c Pficp Hr^pc if ncp pvf rntir^l nfir^n t/~»t Crtil^p tiTYlpc f 
, 111 llllA ^^tlJiC, LlUCh 11 LlhC CALl dUUltlLlUil 1%J1 OlJlI^C lllllCO . 






" fctn it ciTviiilQtp p\?pnt_V\!icpri cti'citpfTipc V 
. Call IL ollllUldLC CVCIIL Ud&CU ollaLCglCa 


yes 




. Ill LIIlo CaoC, io tllC illlCglallUil ocllCIiiC CAcU^l .' 


y Co 


Parflllplism 

L ClL CX.11\jLl^ ILL 


■ HnpQ it *;iinnnrt narnllpl nrnpp'jsiTitr ^ 

• \JIKJ^J 11 i^LILfL/WlL L/LllLlil^l L/l W^L^L^>3lll£i • 


no'' 


Graphics 


! does it have a graphical interface ? 


no, but a programmatic interface 


Simple analysis 


: is it possible to perform simple analysis ? 


yes with visualization 


Complex analysis 


: can more complex analysis be done ? 


itcan^ 


Interface 


: is interface to outside signals possible ? 


indeed'* 




: is it interoperable with other simulators ? 


4 

yes 


Save option 


: can simulation be halted / resumed ? 


yes 


Neuron models 


: can it simulate HH models ? 


itcan^ 




: can it simulate leaky IF models ? 


yes 




: can it simulate multivariate IF models ? 


itcan^ 




: can it simulate conductance-based synaptic interactions ? 


yes 




: can it simulate short-term plasticity ? 


it can^ 




: can it simulate long-term plasticity ? 


itcan^ 




: can it simulate compartmental models with dendrites ? 


no 


Usage 


Development 


: is it currently developed ? 


yes, still a- version 




: how may developers yet ? 


half-time researcher + students 


Support 


: is it supported 


yes 




: what kind of support 


email -I- phone 




: are they user cooperative tools ? 


not yet, tools available® 


Manual 


: are there tutorials and reference material available ? 


yes 




: are there published books on the simulator ? 


no (useless) 




: is there a list of publications of articles that used it ? 


yes 


Import/export 


: is standard (XML) specification import/export available ? 


itcan^ 


Web site 


: is there a web site where all can be found ? 


http : / / enas .gforge.inria.fr 


Source code 


: are there codes available on the web ? 


yes 


Operating system 


: does it run under Linux 


yes, tested 




: does it run under Max-OS X 


yes, tested 




: does it run under Windows 


likely (untested) 


Interoperability 


: using which language can it be used ? 


C/C-I-I-, Python, Java, Php 




: can it be used from other platforms ? 


yes* 



Notes: 

1 : Since clock/event-based mechanisms can be mixed in a event-based simulation, spike times extrapolation is no more to be used, but exact event-time instead. 
2: Exact integration scheme is to be used when allowed by the model, lower-bound spike time evaluation is a new alternative proposed here when the former is not 

possible. 

3: More complex analysis and XML specification import/export is indeed possible, using this kernel within the PyNN environment. The goal was to only develop 
here, what was not available elsewhere. The API has been carefully designed for this puipose. 

4: The interface capability is a key feature of this middle-ware implementation, including in real time applications using spike computations. 

5: Plasticity and other existing models, not discussed here, can be implemented with this middle-ware, and STDP is already considered at that time. The real goal 

is however not to implement "all" models, other simulators do that better, but to propose also alternatives to existing models, as discussed in this paper. 

6: The present development is installed on a forge, thus has all forum/bug-tracking/user-resquest-ticketting, etc.. available. 

7: Though parallelism is not available yet, and the issue not addressed here, the network simulation kernel has been designed and implemented in order to be able 
to cotmect to other kernels, via input/output events. It is thus a feasible extension to run several kernels in parallel, with the drawback that the slower kernel is going 

to drive other kernels local times. 

8: Links with these external platforms such as PyNN (thus NEURON, MvaSpike, ..), NeuroConstruct are made available by the multi-language operability, while 
Scilab and Matlab usage is docimented. 

Table 3: Summary of the main features of the implemented event-based simulation kernel, using the criteria 
proposed to compare existing simulators, see text for details. The required features have been set by the authors 
group of (Brette et aL, 2007) and applied to almost all existing simulators at this date. 



